//*********************************************
//  This is Geant4 Template
//                                  author:Qian
//

#include "MyG4BasedAnalysis.hh"
#include "g4root.hh"
#include "Verbose.hh"
#include <cmath>
#include "G4AntiNeutrinoE.hh"
#include "G4AntiNeutrinoMu.hh"
#include "G4AntiNeutrinoTau.hh"
#include <cstdlib>
#include "G4ProcessType.hh"
#include <ctime>
//fNotDefined, fTransportation, fElectromagnetic, fOptical, fHadronic, fPhotolepton_hadron,
//fDecay, fGeneral, fParameterisation, fUserDefined, fParallel, fPhonon, fUCN
#include "G4DecayProcessType.hh"
//DECAY, DECAY_WithSpin, DECAY_PionMakeSpin, DECAY_Radioactive, DECAY_Unknown, DECAY_External
#include "G4HadronicProcessType.hh"
//fHadronElastic, fHadronInelastic, fCapture, fFission, fHadronAtRest, fLeptonAtRest, fChargeExchange, fRadioactiveDecay
#include "G4TransportationProcessType.hh"
//TRANSPORTATION, COUPLED_TRANSPORTATION, STEP_LIMITER, USER_SPECIAL_CUTS, NEUTRON_KILLER
#include "G4StepStatus.hh"
//fWorldBoundary, fGeomBoundary, fAtRestDoItProc, fAlongStepDoItProc, fPostStepDoItProc, fUserDefinedLimit, fExclusivelyForcedProc, fUndefined
#include "G4TrackStatus.hh"
//fAlive, fStopButAlive, fStopAndKill, fKillTrackAndSecondaries, fSuspend, fPostponeToNextEvent

#include "MyDetectorConstruction.hh"
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wshadow"
#include "TH1F.h"
#pragma GCC diagnostic pop

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

//double MyG4BasedAnalysis::QE[14] = {45., 40., 35., 32., 29., 27., 22., 20., 14., 9., 4., 2., 0.1, 0.0};
MyG4BasedAnalysis::MyG4BasedAnalysis()
{
    SetFileName("g4output.root");

    //-------
    //#ANALYSIS 1. 初始化变量
    fTrkLen = 0;
    depositEnergyTot = 0;
    pTrans = 0;
    fEdeps.clear();
    fHitsX.clear();
    fHitsY.clear();
    fHitsZ.clear();
    pHitsNum = 0;

    motherTrack = NULL;
    //double QE[14]={1,1,1,1,1,1,1,1,1,1,1,1,1,1};
    double QE[25] = {0.34461, 0.428418, 0.3675, 0.2852, 0.2212, 0.1756, 0.1360, 0.1025, 0.07218, 0.05549, 0.04489, 0.03431, 0.02372, 0.01618, 0.008638, 0.0056712, 0.0027077, 0.0028, 0.0013554, 0};
    //double QE[14] = {0.45, 0.40, 0.35, 0.32, 0.29, 0.27, 0.22, 0.20, 0.14, 0.09, 0.04, 0.02, 0.001, 0.0};
    fQE = new TH1F("fQE", "fQE", 25, 120, 210);
    for (int i = 0; i < 25; i++)
        fQE->SetBinContent(i, QE[i]);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

MyG4BasedAnalysis::~MyG4BasedAnalysis()
{
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void MyG4BasedAnalysis::BeginOfRunAction()
{
    if (!active)
        return;

    if (verbose)
        G4cout << "====>MyG4BasedAnalysis::BeginOfRunAction()" << G4endl;

    auto analysisManager = G4AnalysisManager::Instance();
    // Default settings
    analysisManager->SetNtupleMerging(true);
    // Note: merging ntuples is available only with Root output

    analysisManager->SetVerboseLevel(1);
    analysisManager->OpenFile(fFileName);

    //-------
    //#ANALYSIS 2. 定义Ntuple结构

    // Creating 1D histograms
    //analysisManager->SetFirstHistoId(1);
    //analysisManager->CreateH1("phEng", "photon energy", 50, 0., 100); // h1 Id = 0

    // Creating 2D histograms
    //analysisManager->CreateH2("HitOnAnode", "Cherenkov photon hits on the anodes", // h2 Id = 0
    //                          50, -1000., 1000, 50, -1000., 1000.);

    // Creating ntuple
    //
    analysisManager->SetFirstNtupleId(1);
    analysisManager->CreateNtuple("Gun", "Hits");
    analysisManager->CreateNtupleDColumn("X");
    analysisManager->CreateNtupleDColumn("Y");
    analysisManager->CreateNtupleDColumn("Z");
    analysisManager->CreateNtupleDColumn("PX");
    analysisManager->CreateNtupleDColumn("PY");
    analysisManager->CreateNtupleDColumn("PZ");
    analysisManager->CreateNtupleDColumn("kEng");
    analysisManager->CreateNtupleDColumn("PID");
    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("Rec", "Hits");
    analysisManager->CreateNtupleDColumn("Trk1X");
    analysisManager->CreateNtupleDColumn("Trk1Y");
    analysisManager->CreateNtupleDColumn("Trk1Z");

    analysisManager->CreateNtupleDColumn("Trk2X");
    analysisManager->CreateNtupleDColumn("Trk2Y");
    analysisManager->CreateNtupleDColumn("Trk2Z");

    analysisManager->CreateNtupleDColumn("AnodeX");
    analysisManager->CreateNtupleDColumn("AnodeY");
    analysisManager->CreateNtupleDColumn("AnodeZ");

    analysisManager->CreateNtupleDColumn("MirrorX");
    analysisManager->CreateNtupleDColumn("MirrorY");
    analysisManager->CreateNtupleDColumn("MirrorZ");

    analysisManager->CreateNtupleDColumn("Trk3X");
    analysisManager->CreateNtupleDColumn("Trk3Y");
    analysisManager->CreateNtupleDColumn("Trk3Z");

    analysisManager->CreateNtupleDColumn("Trk4X");
    analysisManager->CreateNtupleDColumn("Trk4Y");
    analysisManager->CreateNtupleDColumn("Trk4Z");

    analysisManager->FinishNtuple();

    analysisManager->CreateNtuple("CerenkovPosition", "Hits");
    analysisManager->CreateNtupleDColumn("X", fHitsX);
    analysisManager->CreateNtupleDColumn("Y", fHitsY);
    analysisManager->CreateNtupleDColumn("Z", fHitsZ);
    analysisManager->FinishNtuple();

    return;
}

void MyG4BasedAnalysis::EndOfRunAction()
{
    if (!active)
        return;

    if (verbose)
        G4cout << "====>MyG4BasedAnalysis::EndOfRunAction()" << G4endl;

    //-------
    //#ANALYSIS 6. 在Run结束的时候将ntuple保存到文件
    auto analysisManager = G4AnalysisManager::Instance();
    analysisManager->Write();
    analysisManager->CloseFile();
    return;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void MyG4BasedAnalysis::BeginOfEventAction(const G4Event *)
{
    if (!active)
        return;

    if (verbose)
        G4cout << "====>MyG4BasedAnalysis::BeginOfEventAction()" << G4endl;

    //-------
    //#ANALYSIS 3. 初始化Event开始的参数
    pHitsNum = 0;
    fEdeps.clear();
    fHitsX.clear();
    fHitsY.clear();
    fHitsZ.clear();
    gunPos = G4ThreeVector(.0, .0, .0);
    gunMom = G4ThreeVector(.0, .0, .0);
    fTrack1 = G4ThreeVector(-999, -999, -999);
    fTrack2 = G4ThreeVector(-999, -999, -999);
    fTrack3 = G4ThreeVector(-999, -999, -999);
    fTrack4 = G4ThreeVector(-999, -999, -999);
    fAnode = G4ThreeVector(-999, -999, -999);
    fMirror = G4ThreeVector(-999, -999, -999);

    number = 0;

    motherTrack = NULL;

    return;
}

void MyG4BasedAnalysis::EndOfEventAction(const G4Event *)
{
    if (!active)
        return;

    if (verbose)
        G4cout << "====>MyG4BasedAnalysis::EndOfEventAction()" << G4endl;

    //-------
    //#ANALYSIS 5. 在Event结束的时候将数据保存到ntuple

    auto analysisManager = G4AnalysisManager::Instance();
    analysisManager->FillNtupleDColumn(2, 0, fTrack1.x());
    analysisManager->FillNtupleDColumn(2, 1, fTrack1.y());
    analysisManager->FillNtupleDColumn(2, 2, fTrack1.z());
    analysisManager->FillNtupleDColumn(2, 3, fTrack2.x());
    analysisManager->FillNtupleDColumn(2, 4, fTrack2.y());
    analysisManager->FillNtupleDColumn(2, 5, fTrack2.z());
    analysisManager->FillNtupleDColumn(2, 6, fAnode.x());
    analysisManager->FillNtupleDColumn(2, 7, fAnode.y());
    analysisManager->FillNtupleDColumn(2, 8, fAnode.z());
    analysisManager->FillNtupleDColumn(2, 9, fMirror.x());
    analysisManager->FillNtupleDColumn(2, 10, fMirror.y());
    analysisManager->FillNtupleDColumn(2, 11, fMirror.z());
    analysisManager->FillNtupleDColumn(2, 12, fTrack3.x());
    analysisManager->FillNtupleDColumn(2, 13, fTrack3.y());
    analysisManager->FillNtupleDColumn(2, 14, fTrack3.z());
    analysisManager->FillNtupleDColumn(2, 15, fTrack4.x());
    analysisManager->FillNtupleDColumn(2, 16, fTrack4.y());
    analysisManager->FillNtupleDColumn(2, 17, fTrack4.z());
    analysisManager->AddNtupleRow(2);
    analysisManager->AddNtupleRow(3);
    return;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//此处是产生新的径迹
G4ClassificationOfNewTrack MyG4BasedAnalysis::ClassifyNewTrack(const G4Track *aTrack)
{
    if (!active)
        return fUrgent;

    if (verbose)
        G4cout << "====>MyG4BasedAnalysis::ClassifyNewTrack()" << G4endl;

    //-------
    //#ANALYSIS 4.1 在生成新Track的时候保存相应数据

    G4int parentID = aTrack->GetParentID();

    const G4ParticleDefinition *particle = aTrack->GetParticleDefinition();
    G4int pdgID = particle->GetPDGEncoding();
    G4int charge = particle->GetPDGCharge();
    G4ThreeVector momDir = aTrack->GetMomentumDirection();

    if (charge == 0 && pdgID == 0 && parentID == 1)
    {
        G4String proName = aTrack->GetCreatorProcess()->GetProcessName();
        if (proName == "Cerenkov")
        {
            G4double waveLen = 0.0012398 / aTrack->GetKineticEnergy();
            if (waveLen > 210 || waveLen < 120)
                return fKill;
            double qe = fQE->Interpolate(waveLen);
            if (G4UniformRand() > qe)
            {
                return fKill;
            }
        }
    }

    if (parentID == 0 && aTrack->GetCurrentStepNumber() == 0)
    {
        gunMom = aTrack->GetMomentumDirection();
        gunPos = aTrack->GetPosition();
        G4int pdgID = aTrack->GetParticleDefinition()->GetPDGEncoding();

        auto analysisManager = G4AnalysisManager::Instance();
        analysisManager->FillNtupleDColumn(1, 0, gunPos.x());
        analysisManager->FillNtupleDColumn(1, 1, gunPos.y());
        analysisManager->FillNtupleDColumn(1, 2, gunPos.z());
        analysisManager->FillNtupleDColumn(1, 3, gunMom.x());
        analysisManager->FillNtupleDColumn(1, 4, gunMom.y());
        analysisManager->FillNtupleDColumn(1, 5, gunMom.z());
        analysisManager->FillNtupleDColumn(1, 6, aTrack->GetKineticEnergy());
        analysisManager->FillNtupleDColumn(1, 7, pdgID);
        analysisManager->AddNtupleRow(1);
    }
    return fUrgent;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MyG4BasedAnalysis::PreTrackingAction(const G4Track *aTrack)
{
    if (!active)
        return;

    if (verbose)
        G4cout << "====>MyG4BasedAnalysis::PreTrackingAction()" << G4endl;

    //-------
    //#ANALYSIS 4.2 在Tracking产生的时候保存相应数据

    return;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MyG4BasedAnalysis::PostTrackingAction(const G4Track *aTrack)
{
    if (!active)
        return;

    if (verbose)
        G4cout << "====>MyG4BasedAnalysis::PostTrackingAction()" << G4endl;

    //-------
    //#ANALYSIS 4.3 在Tracking终止的时候保存相应数据
    return;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MyG4BasedAnalysis::SteppingAction(const G4Step *aStep)
{
    if (!active)
        return;

    if (verbose)
        G4cout << "====>MyG4BasedAnalysis::SteppingAction()" << G4endl;

    //---
    //1. 相关参数的获取

    //1.1 Track的相关参数
    G4int parentID = aStep->GetTrack()->GetParentID();
    const G4ParticleDefinition *particle = aStep->GetTrack()->GetParticleDefinition();
    G4int pdgID = particle->GetPDGEncoding();
    G4int charge = particle->GetPDGCharge();
    // G4String particleName = particle->GetParticleName();

    if (parentID == 0 && motherTrack == NULL)
        motherTrack = aStep->GetTrack();

    //以下是G4Track常见的一些参数获取方法
    {
        /*
        G4int iStep = aTrack->GetCurrentStepNumber();
        G4int trackID = aTrack->GetTrackID();
        G4int parentID = aTrack->GetParentID();
        G4ParticleDefinition particle = aTrack->GetParticleDefinition();
        G4int pdgID = particle->GetPDGEncoding();
        G4String pdgName = particle->GetParticleName();
        G4cout << "-->" << iStep << ": tid=" << trackID << " " << pdgName << "(" << pdgID << ") parentID=" << parentID << G4endl;

        G4double energy = aTrack->GetKineticEnergy(); // kinetic energy at current step
        G4double engtot = aTrack->GetTotalEnergy();   // total energy (including m0)
        G4double velocity = aTrack->GetVelocity();    // mm/s
        G4ThreeVector moment = aTrack->GetMomentum();
        G4ThreeVector momDir = aTrack->GetMomentumDirection(); //unit vector
        G4ThreeVector polar = aTrack->GetPolarization();
        G4cout << "-->" << energy << " " << engtot << " " << velocity << " " << moment << " " << momDir << G4endl;

        //two different ways to get the process name
        G4String crtName = (parentID != 0) ? aTrack->GetCreatorModelName() : "none";
        G4String proName = (parentID != 0) ? aTrack->GetCreatorProcess()->GetProcessName() : "none";
        G4cout << "-->" << crtName << " " << proName << G4endl;

        G4ThreeVector trkPos = aTrack->GetPosition(); //current position of this track at the postStepPoint of this step
        G4double trkLen = aTrack->GetTrackLength(); //current track-length of this track at this postStepPoint of this step
        G4double gtime = aTrack->GetGlobalTime();
        G4double ltime = aTrack->GetLocalTime();
        G4double ptime = aTrack->GetProperTime();
        G4cout << "-->trkPos=" << trkPos << " " << trkLen << " " << gtime << " " << ltime << " " << ptime << G4endl;
        */
    }

    //1.2 Step的相关参数
    G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
    G4StepPoint *postStepPoint = aStep->GetPostStepPoint();

    //以下是G4Step常见的一些参数获取方法
    {
        /*
        G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
        G4StepPoint *postStepPoint = aStep->GetPostStepPoint();

        G4double stepLen = aStep->GetStepLength(); //the same as aTrack->GetTrackLength()
        G4double engDep = aStep->GetTotalEnergyDeposit(); // energy deposit inside this step
        G4double engDep2 = aStep->GetNonIonizingEnergyDeposit();
        G4bool isFirst = aStep->IsFirstStepInVolume();
        G4bool isLast = aStep->IsLastStepInVolume();
        G4cout << "-->" << stepLen << " " << engDep << " " << engDep2 << " " << isFirst << " " << isLast << G4endl;

        G4ThreeVector prePos = preStepPoint->GetPosition();
        G4double preGtime = preStepPoint->GetGlobalTime();
        G4double preLtime = preStepPoint->GetLocalTime();
        G4double prePtime = preStepPoint->GetProperTime();
        G4cout << "-->prePos=" << prePos << " " << preGtime << " " << preLtime << " " << prePtime << G4endl;
        G4ThreeVector postPos = postStepPoint->GetPosition(); // same as the information get from aTrack
        G4double postGtime = postStepPoint->GetGlobalTime();
        G4double postLtime = postStepPoint->GetLocalTime();
        G4double postPtime = postStepPoint->GetProperTime();
        G4cout << "-->postPos=" << postPos << " " << postGtime << " " << postLtime << " " << postPtime << G4endl;
        
        G4ThreeVector preMomDir = preStepPoint->GetMomentumDirection();
        G4ThreeVector preMomtum = preStepPoint->GetMomentum();
        G4double preEngTot = preStepPoint->GetTotalEnergy();
        G4double preEng = preStepPoint->GetKineticEnergy();
        G4double preV = preStepPoint->GetVelocity();
        G4double preBeta = preStepPoint->GetBeta();
        G4double preGamma = preStepPoint->GetGamma();

        G4ThreeVector postMomDir = postStepPoint->GetMomentumDirection();
        G4ThreeVector postMomtum = postStepPoint->GetMomentum();
        G4double postEngTot = postStepPoint->GetTotalEnergy();
        G4double postEng = postStepPoint->GetKineticEnergy();
        G4double postV = postStepPoint->GetVelocity();
        G4double postBeta = postStepPoint->GetBeta();
        G4double postGamma = postStepPoint->GetGamma();
        
        G4double mass = preStepPoint->GetMass();
        G4double charge = preStepPoint->GetCharge();
        G4double magmom = preStepPoint->GetMagneticMoment();

        auto *pVolume = postStepPoint->GetTouchableHandle()->GetVolume();
        auto *pLogicVolume = (pVolume!=NULL) ? pVolume->GetLogicalVolume() : NULL;

        */
    }

    //1.3 以下是拿DetectorConstruction相关参数的方法，根据具体情况修改即可
    {
        /*
        //拿Detector指针
        const MyDetectorConstruction* detectorConstruction = static_cast<const MyDetectorConstruction*> (G4RunManager::GetRunManager()->GetUserDetectorConstruction());

        //查找已知名字的逻辑体
        G4Box *fBox; 
        G4LogicalVolume* fLogicVol = G4LogicalVolumeStore::GetInstance()->GetVolume("YourLogicVolumeName");
        if ( fLogicVol ) fBox = dynamic_cast<G4Box*>(fLogicVol->GetSolid());
        fBox->GetXHalfLength() //这样就可以拿到具体的参数了。根据具体情况修改即可
        */
    }

    //-------
    //#ANALYSIS 4.4 在Steppinging的时候保存相应数据

    //---
    //2. 添加一些判断，并保存对应的数据。以下为演示，且按ANALYSIS 2. Ntuple定义的结构进行保存

    //G4String proName = postStepPoint->GetProcessDefinedStep()->GetProcessName();
    //G4ProcessType proType = postStepPoint->GetProcessDefinedStep()->GetProcessType();
    //G4int proSubType = postStepPoint->GetProcessDefinedStep()->GetProcessSubType();
    auto *preVolume = preStepPoint->GetTouchableHandle()->GetVolume();
    auto *postVolume = postStepPoint->GetTouchableHandle()->GetVolume();
    G4String PreVolume = preVolume->GetName();
    //为了获取前、后步所在的volume的名称
    if (postVolume == NULL)
        return;
    G4String name = particle->GetParticleName();
    G4String proName = postStepPoint->GetProcessDefinedStep()->GetProcessName();
    //G4int parentID = aTrack->GetParentID();
    G4LogicalVolume *preNowVolume = preVolume->GetLogicalVolume();
    G4LogicalVolume *postNowVolume = postVolume->GetLogicalVolume();
    G4ThreeVector postPos = postStepPoint->GetPosition();
    G4double optX = postPos.x();
    G4double optY = postPos.y();
    G4double optZ = postPos.z();
    G4ThreeVector postMom = postStepPoint->GetMomentum();
    G4double momX = postMom.x();
    G4double momY = postMom.y();
    G4double momZ = postMom.z();
    G4double momTot = sqrt(momX * momX + momY * momY + momZ * momZ);
    G4ThreeVector postMomDir = postStepPoint->GetMomentumDirection();
    G4ThreeVector preMomDir = preStepPoint->GetMomentumDirection();
    G4double preEng = preStepPoint->GetKineticEnergy();
    G4double postEng = postStepPoint->GetKineticEnergy();
    //G4double waveLen = 0.0012398 / aTrack->GetKineticEnergy();
    G4double waveLen = 0.001242375 / postEng;

    auto analysisManager = G4AnalysisManager::Instance();

    G4String particleName = aStep->GetTrack()->GetDynamicParticle()->GetDefinition()->GetParticleName();

    if (particleName == "opticalphoton") // for the 'opticalphoton'
    {
        if (aStep->GetTrack()->GetTrackStatus() != fAlive)
        {

            G4ThreeVector start = aStep->GetTrack()->GetVertexPosition();
            G4ThreeVector momentum = aStep->GetTrack()->GetVertexMomentumDirection();
            G4ThreeVector final = aStep->GetTrack()->GetPosition();

            //if (_z < -999 && _z > -1001)
            {
                fHitsX.push_back(final.x());
                fHitsY.push_back(final.y());
                fHitsZ.push_back(final.z());
                //analysisManager->FillNtupleDColumn(2, 0, final.x());
                //analysisManager->FillNtupleDColumn(2, 1, final.y());
                //analysisManager->FillNtupleDColumn(2, 2, final.z());
                //analysisManager->AddNtupleRow(1);
            }
            return;
        }
    }

    if (postNowVolume->GetName() == "track1Vol" && parentID == 0 && fTrack1.x() == -999 && fTrack1.y() == -999 && fTrack1.z() == -999)
    {
        fTrack1 = aStep->GetTrack()->GetPosition();
    }

    if (postNowVolume->GetName() == "track2Vol" && parentID == 0 && fTrack2.x() == -999 && fTrack2.y() == -999 && fTrack2.z() == -999)
    {
        fTrack2 = aStep->GetTrack()->GetPosition();
    }

    if (postNowVolume->GetName() == "track3Vol" && parentID == 0 && fTrack3.x() == -999 && fTrack3.y() == -999 && fTrack3.z() == -999)
    {
        fTrack3 = aStep->GetTrack()->GetPosition();
    }

    if (postNowVolume->GetName() == "track4Vol" && parentID == 0 && fTrack4.x() == -999 && fTrack4.y() == -999 && fTrack4.z() == -999)
    {
        fTrack4 = aStep->GetTrack()->GetPosition();
    }

    if (postNowVolume->GetName() == "boardVol" && parentID == 0 && fAnode.x() == -999 && fAnode.y() == -999 && fAnode.z() == -999)
    {
        fAnode = aStep->GetTrack()->GetPosition();
    }

    if (postNowVolume->GetName() == "ballVol" && parentID == 0 && fMirror.x() == -999 && fMirror.y() == -999 && fMirror.z() == -999)
    {
        fMirror = aStep->GetTrack()->GetPosition();
    }
}
